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""^3 ■ Abstract 

Qh The spherical tensor gradient operator y™(V), which is obtained by re- 

fH , placing the Cartesian components of r by the Cartesian components of 

V in the regular solid harmonic yf 1 (r), is an irreducible spherical tensor 

S of rank £, Accordingly, its application to a scalar function produces an 
irreducible spherical tensor of rank I. Thus, it is in principle sufficient 
;> to consider only multicenter integrals of scalar functions: Higher angu- 

lar momentum states can be generated by differentiation with respect 
to the nuclear coordinates. Many of the properties of y™(V) can be 
understood easily with the help of an old theorem on differentiation by 
Hobson [Proc. London Math. Soc. 24, 54 - 67 (1892)]. It follows from 
Hobson's theorem that some scalar functions of considerable relevance 
as for example the Coulomb potential, Gaussian functions, or reduced 
Bessel functions produce particularly compact results if y™(V) is ap- 
plied to them. Fourier transformation is very helpful to understand 
the properties of y" l iy) since it produces y™ (— ip). It is also possi- 
ble to apply yj"(V) to generalized functions, yielding for instance the 
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spherical delta function 5|™(r). The differential operator y™(V) can 
also be used for the derivation of pointwise convergent addition theo- 
rems. The feasibility of this approach is demonstrated by deriving the 
addition theorem of r v y™ir) with u G R. 
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1 Introduction 

In the early 17th century, scientific methodology advanced significantly by 
what Kline |561 Chapter 16] called the mathematization of science. Roughly 
at the same time, the foundations for the later development of calculus were 
laid (see for example |56l Chapter 17]). Differential and integral calculus 
greatly extended the arsenal of mathematical techniques that could be used 
for a description and analysis of scientific phenomena. Ultimately, this de- 
velopment led to a period of unity between mathematics and sciences which 
lasted approximately to the end of the 19th or to the beginning of the 20th 
century. In that period, it was frequently not possible to decide whether 
somebody was predominantly a mathematician or predominantly a scien- 
tist, and mathematics and the sciences enriched each other greatly by cross 
fertilization. 

A striking example for this unity between mathematics and the sciences 
is provided by Peter Debye who had been a student of Sommerfeld. Among 
chemists, Debye is best known for his work on electrolyte solutions, which 
earned him the Nobel Prize in Chemistry in 1936. It is, however, not nearly 
so well known that he was also an excellent mathematician who made sig- 
nificant contribution to the theory of Bessel and Hankel functions (see for 
example |115( Chapter VIII]). 

Since the times of Debye, the amount of collective knowledge has in- 
creased tremendously. Accordingly, contemporary research is predominantly 
done by specialists, who know almost everything about almost nothing, 
whereas generalists like Debye, who had done excellent research in math- 
ematics, physics, and physical chemistry, have become exceedingly rare. In 
particular, there is a widening gap between mathematicians and those that 
apply mathematics. In the past, mathematicians had developed new ana- 
lytical or numerical techniques that greatly helped scientists and engineers 
to do research in their disciplines. In return, open scientific or engineering 
problems had always provided a valuable source of inspiration for mathe- 
maticians. Unfortunately, this is now longer true. Due to increased special- 
ization, communication between mathematics and those disciplines, that use 
mathematics as their main language, has deteriorated. This is certainly a 
deplorable development because cross fertilizations typically happens at the 
interfaces of different disciplines. 

Fortunately, encouraging counterexamples still exist: For example, Pro- 
fessor Josef Paldus started his career at the Heyrovsky Institute in Prague 
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doing experimental and theoretical work in electrochemistry, and he later 
also worked as an experimental spectroscopist. In addition, he soon ven- 
tured into the newly emerging field of quantum chemistry. In 1968, he came 
as Brezhnev's gift to the Department of Applied Mathematics of the Uni- 
versity of Waterloo. There, his predominant research interest has been the 
treatment of electronic correlation, emphasizing highly sophisticated math- 
ematical techniques as for example the representation theory of Lie groups. 
It was this strong mathematical orientation of his research which earned him 
- in addition to numerous other honors - the Fellowship of the prestigious 
Fields Institute for Research in Mathematical Sciences. 

Quite in the spirit of the highly interdisciplinary research of Professor 
Josef Paldus, I want to discuss in this article a topic that is located some- 
where at the interface between mathematics and atomic and molecular elec- 
tronic structure theory. 

Since space is isotropic, free atoms are spherically symmetric and an- 
gular momentum is conserved. Thus, spherical polar coordinates and the 
machinery of angular momentum theory lead to significant computational 
and conceptual simplifications in atomic structure calculations. It is prob- 
ably only a slight exaggeration to claim that atomic electronic structure 
theory is essentially an application of angular momentum theory. 

In the case of molecules, the benefits of spherical coordinates and angular 
momentum theory are not so obvious. In general, molecules possess no 
spatial symmetry at all, and if they have one, it is a lower symmetry than 
spherical symmetry. Thus, neither spherical polar coordinates nor angular 
momentum theory lead to a such spectacular reduction of computational 
complexity as they do in the case of atoms. 

However, chemists have found it helpful to approach molecular electronic 
structure theory from the perspective of atoms. For example, the most suc- 
cessful computational scheme, the so-called LCAO-MO approach, is based 
on the tacit assumption that the parentage of atoms facilitates our attempts 
of understanding the electronic structure of molecules. Thus, in molecular 
electronic structure calculations it is common to use angular momentum 
theory at least locally, i.e., with respect to the atomic centers. In this way, 
at least some of the formal advantages of angular momentum theory can be 
retained. 

Nevertheless, the use of angular momentum theory for systems lacking ro- 
tational symmetry causes new mathematical and computational challenges. 
In the case of atoms, only relative small angular momentum quantum num- 
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bers occur, and the coupling of irreducible spherical tensors produces finite 
expressions consisting of a few terms only. In the case of molecules, the 
situation is in general much more complicated: We have to deal with infi- 
nite series expansions over angular momentum states that do not necessarily 
converge rapidly. Thus, we must be able to compute the typical quantities 
of angular momentum theory both efficiently and reliably even for very large 
angular momentum quantum numbers. Quite often, this is not so easy. 

In molecular electronic structure calculations with analytic basis func- 
tions centered at the atomic nuclei, the basic quantities are matrix elements, 
i.e., essentially multicenter integrals whose evaluation may be very difficult. 
It is an undeniable empirical fact that it is often comparatively easy to ob- 
tain explicit analytical expressions for multicenter integrals over functions 
that are scalars or irreducible spherical tensors of rank zero with respect to 
their local (atomic) coordinate systems. If, however, the functions occurring 
in the multicenter integral are irreducible spherical tensors of higher ranks, 
one can easily get lost in an algebraic jungle and the derivation of explicit 
expressions can become extremely difficult or even impossible. 

It is another empirical fact that it is usually much easier to differentiate 
than to integrate. Accordingly, it is an obvious idea to try to generate an 
explicit expression for a multicenter integral over nonscalar functions by dif- 
ferentiating the simpler expression for the corresponding integral over scalar 
functions (preferably the simplest scalar functions) with respect to scaling 
parameters and/or nuclear coordinates (compare also [1221 Section IV]). The 
use of generating differential operators does not necessarily produce closed 
form expression that hold for arbitrary quantum numbers. Instead, it may 
be necessary to derive a new expression for each set of quantum numbers. 
Obviously, a multitude of special formulas is not nearly as convenient as a 
neat general explicit expression, but powerful computer algebra systems like 
Maple or Mathematica with the ability of automatically generating FOR- 
TRAN or C code help to make such an approach feasible (see for example 

USD- 

It is relatively easy to generate multicenter integrals of higher scalar func- 
tions by differentiating the most simple scalar functions with respect to their 
scaling parameters. In the case of a Is Slater-type or Gaussian function, we 
can construct higher scalar functions easily by repeatedly using the relation- 
ships <9exp(— ar)/da = — rexp(— ar) or <9exp(— ar 2 )/da = — r 2 exp(— ar). 

The generation of anisotropic functions, that are irreducible spherical 
tensors of rank £, from scalar functions is less straightforward: It follows 
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from Hobson's theorem |42j . which is discussed in Section |3J that we have 
to apply the spherical tensor gradient operator y™{V) to scalar functions. 
This differential operator is obtained by replacing in the regular solid har- 
monic y™{r) = r YF^iQ, 4>) the Cartesian components of r = (x, y, z) by the 
Cartesian components of V = (d/dx,d/dy,d/dz). 

I came across the differential operator y™(V) during my PhD thesis |117j 
in which I worked on multicenter integrals of reduced Bessel functions and 
their anisotropic generalizations, the so-called B functions. These functions, 
which are discussed in Sectional are a special class of exponentially decaying 
functions with some very useful properties. When I derived the remarkably 
compact Fourier transform (|5.6|) of a B function |117[ Eq. (7.1-6) on p. 160], 
I noticed that the Fourier integral representation (|5.<S|) of an anisotropic B 
function B™ e (a,r) differs from that of the scalar B function 5^_^ (a,r) 
only by an additional regular solid harmonic y™(— ip)- Because of (|4.4|) . 
3^™(ip) can De interpreted to be the Fourier transform of the differential 
operator ^"(V), and I immediately deduced that the anisotropic B function 
B™ £ (a,r) can according to (|5.9|) be generated by applying y" l {S7) to the 
scalar B function B™ +efi (a,r) [HZ1 Eq. (7.1-10) on p. 161]. 

Of course, this observation aroused my interest: I wanted to know whether 
this result can also be proved directly without the help of Fourier transfor- 
mation, and I also wanted to know whether the differential operator ^^(V) 
is a useful analytical tool in different contexts. After the completion of my 
PhD thesis, I studied the spherical tensor gradient operator more seriously. 
I soon learned that I was not the only one and in particular not the first 
one who had studied and used this differential operator: The first article 
dealing with the spherical tensor gradient operator, which I am aware of, 
was published by Hobson [32] in 1892. I also noticed that y™(V) is a highly 
useful analytical tool for a wide range of problems. For example, I found 
recent articles which describe successful applications of y^iy) in scattering 
[Tl| or in solid state theory |HH] • 

From the perspective of quantum chemistry, it is probably more inter- 
esting that y™(V) can be extremely useful in the context of molecular mul- 
ticenter integrals of exponentially decaying functions, as shown in articles 
by Grotendorst and Steinborn |39| 140] . Niukkanen |70 ( I71 ( 173]. Novosadov 
[13 EE 112 EB1 EE EH, Tai PH], and in my own research [11111121111122111201 

\rmrnm\rm\m\. 

The spherical tensor gradient operator y™(V) is - as discussed in more 
details in Section [31 - also extremely useful in the context of multicenter in- 
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tegrals of spherical Gaussian functions. For example, Dunlap reformulated 
in a series of recent articles [221 1231 1241 1251 126j the theory of multicenter in- 
tegrals of anisotropic spherical Gaussians by systematically applying 3^™(V) 
to the corresponding integrals of Is Gaussians. Many more references on 
multicenter integrals of spherical Gaussians can be found in Section [31 

These examples should suffice to show that yjpCV) is indeed a very useful 
analytical tool that has been applied successfully to a wide range of prob- 
lems. Nevertheless, a reasonably comprehensive discussion the mathematical 
properties of the spherical tensor gradient operator from the perspective of 
a potential future user seems to be missing. This is what I intend to provide 
with this article. 

In Section the spherical tensor gradient operator is introduced and 
some general features - in particular those based on its tensorial nature - 
are discussed. 

Section |31 treats Hobson's theorem on differentiation, by means of which 
many properties of 3^™(V) can be derived and understood easily. 

Fourier transformation, which is discussed in Section 0J maps the gradi- 
ent operator V to ip, or equivalently y™{S7) to y™(yp). Thus, in momen- 
tum space, we only have to study the comparatively simple multiplicative 
operator y^iip) whose properties are very well understood. This greatly 
facilitates a theoretical analysis. 

In Section [51 the so-called reduced Bessel functions and their anisotropic 
generalizations, the so-called B functions, are discussed. These functions 
play a special role among exponentially decaying functions because of their 
exceptionally simple Fourier transform and because of the ease, with which 
the spherical tensor gradient operator can be applied. 

Classically, the domain of the spherical tensor gradient operator consists 
of the differentiable functions / : M 3 — > C, but it makes sense to apply it 
also to generalized functions. Thus, in Section H3 the spherical delta function 
and related objects - for example distributional B functions - are treated. 

In Section [71 the derivation of addition theorems of essentially arbitrary 
irreducible spherical tensors with the help of the spherical tensor gradient 
operator is discussed. In addition, the addition theorem of v^y^ir) with 
v € R is constructed in order to demonstrate the feasibility of the whole 
approach. 

This article is concluded by a short summary in Sectional Finally, there 
are three Appendices: The terminology used in this article is introduced in 
Appendix ^ the for our purposes most relevant properties of the spherical 
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harmonics are reviewed in Appendix [Bl and Gaunt coefficients are discussed 
in Appendix O 



2 The Spherical Tensor Gradient Operator 

The spherical tensor gradient operator yjpfV) can be introduced via the 
explicit expression (jB.121) for the regular solid harmonic 3^™(r). This explicit 
expression holds for the Cartesian components of essentially arbitrary three- 
dimensional vectors. Thus, we obtain the differential operator y™(V) if we 
replace in ()B.8|) the Cartesian components of r = (x, y, z) by the Cartesian 
components of V = (d/dx, d/dy, d/dz): 



3T(V) 



21 + 1 
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(2.1) 



A new differential operators is not necessarily a useful thing, let alone a 
major achievement. In atomic and molecular calculations, we are interested 
in functions that are in defined in terms of spherical polar coordinates and 
that can be expressed as a radial part multiplied by a spherical harmonic: 



*T(r) = h(r)Yr(r/r) 



(2.2) 



The straightforward differentiations of such an irreducible tensor of rank I 
with respect to the Cartesian components x, y, and z of r would in general 
lead to extremely messy expressions. Because of the convenient orthonor- 
mality properties of the spherical harmonics it would be advantageous if the 
angular part of such a product could be expressed in terms of spherical har- 
monics. Of course, it should be possible to do the necessary algebra, but 
in particular for large angular momentum quantum numbers we would be 
confronted with nontrivial technical problems. 

Thus, we arrive at the paradoxical statement that the differential opera- 
tor ypCV) is practically useful only if it is not necessary to do differentiations 
with respect to x, y, and z via the defining explicit expression ()2.1|) . For- 
tunately, these differentiations can be avoided since 3^™(V) is just like the 
corresponding regular solid harmonic 3^™(r) an irreducible spherical tensor 
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of rank I (compare j!21 p. 312]). Consequently, matrix elements involving 
3^™(V) and other irreducible spherical tensors can be handled via the pow- 
erful machinery of angular momentum coupling. 

As discussed in more details later, the product y^ 1 (V)F^ 12 (r) can be 
expressed as a finite linear combination of Gaunt coefficients defined in ()C.1|) , 
radial functions 'fy^Jj), and spherical harmonics |1311 Eq. (4.7)]: 

^(V) *£»(r) 

p 

'-max 

= Yl (2) + m 2 \£ imi \e 2 m 2 ) 7 | ll2 (r) Y™ 1+m2 (r/r) . (2.3) 

The summation limits £ m - m and £ ma , x are given by (|C.5|) . 

It is possible to derive alternatives to ()2.3j) which also take into account 
the tensorial nature of both 3^™ 1 (V) and F^ 2 (r). For example, Bayman 
[0] derived an equivalent expression involving Clebsch-Gordan coefficients 
instead of Gaunt integrals as in (|2.3|) . However, it was argued in [1311 
Section III] that V behaves like r under reflections. Accordingly, parity 
is conserved, and the replacement of Clebsch-Gordan coefficients by Gaunt 
integrals according to (|C.2|) leads to a considerable formal simplification. 

As discussed in more details in Section |IJ the functions t| ^ (r) in (|2.3|) 
can be obtained by differentiating the radial part fe 2 ( r ) of the spherical 
tensor F^ 2 (r) with respect to r (compare the discussion in |1311 Sections 
III and IV]). Thus, the systematic exploitation of the tensorial nature of 
y™ 1 (V) makes it possible to replace potentially troublesome differentiations 
with respect to the three Cartesian components x, y, and z by comparatively 
benign differentiations with respect to the radial variable r. This is a very 
important technical aspect since it greatly facilitates the application of the 
spherical tensor gradient operator to spherical tensors of the type of (|2.2jl . 

Other expressions for the product y 1 ^ 1 (V)F^ 2 (r) can be found in articles 
by Santos |25|, Stuart [m], Niukkanen' |23 , and Rashid [HS1- 

3 Hobson's Theorem on Differentiation 

The first article on the spherical tensor gradient operator y^CV), which 
I am aware of, is due to Hobson |42| p. 67] who derived in 1892 a very 
consequential theorem on the differentiation of functions / : W 1 — > C. This 
theorem is also discussed in Hobson's book |43| pp. 124 - 129] that was first 
published in 1931. 
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In contemporary articles dealing with the spherical tensor gradient op- 
erator, Hobson's theorem - Eq. Q3.1j) below - is in most cases completely 
ignored, but in the article by Bott, Methfessel, Krabs, and Schmidt |13l Sec- 
tion IV], the simplified form (j3.2|) of Hobson's theorem was even rederived 
in a less general form. This is an undeserved neglect. Firstly, many prop- 
erties of the spherical tensor gradient operator 3^™(V) can be deduced and 
understood easily via Hobson's theorem. Secondly, it should be possible to 
generalize the approach described in this article, which is based on the three- 
dimensional regular solid harmonics y™{r), to differential operators based 
on n-dimensional hyperspherical harmonics. A problem of obvious physical 
relevance would be the construction of differential operators associated to 
the four-dimensional hyperspherical harmonics, which are for instance dis- 
cussed in books by Avery 013] and Judd [23 Sections 2 and 3 and Appendix 
2] or also in [TTE1 Section VI]. 

In three-dimensional form - which is all we need here - Hobson's theorem 
can be formulated as follows (compare also [T29, Eq. (4.13)]): 
Let f n (x,y,z) be a homogeneous polynomial of degree n G N in the vari- 
ables x,y,z, and let F be any (differentiable) function that depends only 
on r 2 = x 2 + y 2 + z 2 . Then, the application of the differential operator 
fn(d/dx, d/dy, d/dz) to F can be expressed in closed form according to 



Let us now assume that f n is not only a homogeneous polynomial of degree n, 
but also a solution of the three-dimensional homogeneous Laplace equation 
(|B.8|) . Then, (|3.1j) simplifies considerably because in the sum on the right- 
hand side only the power V 2t/ with v = produces a nonzero result: 



As discussed in Appendix^ the regular solid harmonic 3^™( r ) i s a P or y~ 
nomial solution of the homogeneous Laplace equation (jB.8|) of degree i. 
Accordingly, Hobson's theorem applies in its simplified form (|3.2j) . and we 
obtain (compare |129| Eq. (3.11)]): 



Ernst Joachim Weniger: The Spherical Tensor Gradient Operator 




(3.1) 




(3.2) 




(3.3) 



3: Hobson's Theorem on Differentiation 



11 



If we now set F(r 2 ) = tp(r) and use d/dr 2 = l/(2r) d/dr, we obtain (compare 
[1231 Eq. (3.12)]): 

y?{r). (3.4) 

It follows from either (|3.3|) or Q3.4|) that 3^(V) can be viewed as a generat- 
ing differential operator which transforms a scalar function - an irreducible 
spherical tensor of rank - to an irreducible spherical tensor of rank I. 

Next, let us apply the spherical tensor gradient operator y™ 1 (V) to 
irreducible spherical tensors Fj^ 2 of the type of (|2.2|) with nonzero rank l 2 . 
To simplify things, let us for the moment assume that such a spherical tensor 
Fj^ 2 (r) satisfies a relation of the kind of (|3.4j) . i.e., it can be generated by 
applying 3^™ 2 (V) to a suitable scalar function $^ 2 (r): 

F^ 2 (r) = y% 2 (V)* h (r). (3.5) 

In view of y^(V)F^ 2 (r) = y 1 ^ 1 (V)y™ 2 (V) $* a (r), we need an explicit 
expression of manageable complexity for the product y™ 1 (V)3^™ 2 (V). This 
can be accomplished easily. If we multiply the linearization formula l)C.4|) of 
the surface harmonics by r^ 1+ ^ 2 , we obtain the linearization formula of the 
regular solid harmonics: 

^max 

= (2) ( im i + m 2 |4mi |£ 2 m 2 ) r 2A£ ^+ m 2 ( r ) . ( 3 . 6 ) 

It follows at once from the summation limits (fC>.f>|> that the abbreviation A£ 
defined by ()C.6|) is either a positive integer or zero. Accordingly, Q3.6|) can 
also be interpreted as a linearization formula for certain polynomials in the 
Cartesian components of an essentially arbitrary three-dimensional vector r. 
Thus, ()3.6|) also holds for r = V, and we obtain (see for example [701 Eq. 
(15)], [ZS1 Eq. (65)], or [Ml Eq. (4.24)]): 

^?(v)^ 2 ( v ) 

^max 

= (2) (^i+m 2 \£ 1 m 1 \£ 2 m 2 )V 2M y^ +m2 (V). (3.7) 



IdV ( , 
r dr J 
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If we now combine (|3.5|) and (|3.7[) , we obtain an expression which is obviously 
of the form of (|2.3|) (see for example [1331 Eq. (3.9)]): 

^max 

y™ 1 (V) (r) = ^ + m 2 \hmi \hm 2 ) 

P—P ■ 

1 d 



x V 2Ae 



r dr ' 2W 



^ mi+m2 (r). (3.8) 



There are some radially symmetric functions of considerable physical rel- 
evance which produce results of remarkable simplicity if 3^™(V) is applied 
to them via (|3.4|h The classic example is the Coulomb potential 1/r. Hob- 
son |42j showed that the irregular solid harmonic Z™(r) defined in 1)13,71) is 
generated by applying 3^™(V) to 1/r (further details can be found in Hob- 
son's book [1*3 pp. 124 - 129]). In modern notation, Hobson's result can be 
expressed as follows (see for example |129[ Eq. (4.16)]): 

It is possible to derive this relationships differently (see for example |12| 
Chapter 6.18]), but its derivation via ()3.4j) is in my opinion much more 
straightforward and transparent. 

It is also quite easy to obtain an explicit expression for the product 
J™ 1 (V)2™ 2 (r). With the help of (ESJ) and we immediately obtain 

the following compact expression involving Gaunt coefficients, Pochhammer 
symbols, and irregular solid harmonics [1331 Eq. (4.7)]: 

y^(X7)Z™*(r) = (e 1 +£ 2 m 1 + m2 \£ 1 m 1 \e 2 m 2 ) 

x (_2)* (y^i+fr Z„ m i+ m2 (r) . (3.10) 

Alternative derivations of this result are discussed in |12U1 Section 7]. 

In the vast majority of all LCAO-MO electronic structure calculations, 
Gaussian functions are used as basis functions. Traditionally, Cartesian 
Gaussian functions x l y^ z k exp(— ar 2 ) with k € No and a € R+ have dom- 
inated, which indicates that the spherical tensor gradient operator should 
be irrelevant for calculations with Gaussian basis functions. However, in 
recent years research on molecular integrals of Gaussians has increasingly 
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emphasized Gaussian functions whose angular parts are spherical harmon- 
ics. For these functions, 3^™(V) is again an extremely useful mathemati- 
cal tool (see for example the articles by Arakane and Matsuoka pQ, Chow 
Chiu and Moharerrzadeh [El EH CHI, Dimlap [22] E3 EH E3 EI], Fieck 
[2*71 l2l"] . Fortunelli and Carrravetta [32], Fortunelli and Salvetti [?""""], Fu- 
jimura and Matsuoka 34 , Hu, Staufer, Birkenheuer, Igoshine, and Rosch 
[UJ, Ishida HH1 SHI inOl EU E21, Kuang and Lin [SHIES], Maretis """2], Mat- 
suoka [S3 El E23 EHl E7J , and Saunders [H] ) . 

The usefulness of 3^ n (V) in connection with Gaussian functions follows 
at once from the fact that exp(— ar 2 ) is an eigenfunction of the differential 
operator (l/r)d/dr with eigenvalue —2a. Thus, the application of 3^"(V) to 
exp(— ar 2 ) via ()3.4)) yields according to Fieck [2H] the following remarkably 
compact irreducible spherical tensor of rank t: 

3> € m (V) exp(-ar 2 ) = (-2a)' exp(-ar 2 ) yf (r) . (3.11) 

A detailed discussion of the usefulness of the differential operator 3^™(V) 
in the context of Gaussian functions and their multicenter integrals would 
certainly be of considerable interest, but unfortunately it would clearly be 
beyond the scope of this article. In addition, I cannot claim that molecular 
integrals of Gaussian functions are my field of expertise. Therefore, the 
interested reader is referred to the articles listed above. 

Another class of functions, which have close ties to the differential opera- 
tor 3^(V) 3 are the s-called B functions which will be discussed in a detailed 
way in Section [""] 

4 Fourier Transformation 

As discussed in Section Hobson's theorem implies that a spherical tensor 
of rank t can be generated by applying 3^ n (V) to a spherically symmetric 
function <p(r) according to (|3.4|) . If an irreducible spherical tensor F^ 2 of 
the type of (|2.2|) also satisfies (|3.5j) . then the simple explicit expression (|3.8|) 
for the product y™ 1 (V)F^ 2 (r) can be derived via (""""!)) . 

In principle, (|3.8|) should suffice for our purposes, since the scalar function 
<3?^ 2 (r) in 1)3.5)) can according to 1)3.4)1 be obtained from the scalar function 
fe 2 ( r ) m Q2.2)) by repeated integration with respect to r. However, repeated 
integrations are at least potentially a source of trouble, and it is thus desir- 
able to express the product y™ 1 (V)-F^ 12 (r) according to ()2.3)) in terms of 
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radial functions £ (r) that can be obtained by differentiating the radial 
function fg 2 (r) in (|2.2[) . That this is indeed possible can be shown most 
easily with the help of Fourier transformation, which will be used in its sym- 
metrical form. Thus, a function / : M 3 — » C and its Fourier transform / are 
connected by the integrals (|A.2|) and (|A.3|) . 

The practical usefulness of Fourier transformation is obvious, and it 
would clearly be beyond the scope of this article to mention all success- 
ful scientific applications. Let me just mention that Fourier transformation 
is - as first shown by Prosser and Blanchard jH2j and by Geller - one of 
the principal methods of handling molecular multicenter integrals. 

For example, the convolution of two functions /, g : M 3 -> C can be 
expressed as an inverse Fourier integral (see for example j 1 3 H Eq. (5.11)]): 

J f(r-r')g(r')d 3 r' = J f(p) g(p) d 3 p , (4.1) 

We immediately obtain this result if we replace in (|4.1|) / and g by their 
inverse Fourier integrals according to (|A,3|) and use the following integral 
representation of the three-dimensional delta function (see for example [551 
Eq. (3.16)]): 

<5(r) = (2^)- 3 J e ir - p d 3 p. (4.2) 

Six-dimensional integrals describing the Coulomb interaction of two charge 
densities [/(fi)]* and g(r2) can also be expressed as three-dimensional in- 
verse Fourier integrals (see for example [126, Eq. (2.22)]): 

/ / ^ r ^* \ ri -l 2 -R\ ^)dVd 3 r 2 

= 4vr J ^-^[/-(p)]*5(p)d 3 p, (4.3) 

For a proof of (|4.3|) , we also need the Fourier transform 1)6.4(1 of the Coulomb 
potential, which - as discussed in more details in Section 0- holds only in 
the sense of distributions. 

In connection with multicenter integrals in general and with the spheri- 
cal tensor gradient operator in special, Fourier transformation suffers from 
a serious limitation which must not be ignored. Classically, Fourier trans- 
formation is defined only for functions that are absolutely integrable, i.e., 
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which belong to the function space L 1 (M 3 ), but by means of a suitable lim- 
iting procedure it can be extended uniquely to give a unitary map from the 
Hilbert space L 2 (R 3 ) of square integrable functions onto itself |84l Theorem 
IX. 6 on p. 10]. Unfortunately, not all functions of interest are square inte- 
grable and thus possess Fourier transforms that are meaningful in the sense 
of classical analysis. As discussed in more details in Section H3 the Coulomb 
potential and the irregular solid harmonic, which are connected via (|3.9[) . 
possess Fourier transforms only in the sense of generalized functions or dis- 
tributions. For the sake of conceptual simplicity, let us tacitly assume that 
all Fourier integrals in this Section are meaningful in the sense of classical 
analysis. 

Fourier transformation greatly simplifies the treatment of the spherical 
tensor gradient operator because of the obvious relationship 

^ m (V)e irp = i e yp(p)e ir - p . (4.4) 

Thus, in momentum space 3^™(V) is a simple multiplicative operator with 
known coupling properties. Moreover, (|4.4j) makes it plausible that y™(S7) 
is indeed an irreducible spherical tensor of rank £ (compare |12l p. 312]). 

If we use in the Fourier integral (|A.2I) both ()2.2|) as well as the Rayleigh 
expansion of a plane wave in terms of spherical Bessel functions and spherical 
harmonics (compare for instance |12l p. 442]), 

oo £ 

e ±i P r = 4n J2(±l) l Mpr) £ [Yr( P /p)]*Yr(r/r), (4.5) 
1=0 m=-£ 

we see that the Fourier transform F^ 2 (p) of the irreducible spherical tensor 
F^ 2 (r) is also a spherical tensor of rank £ 2 in momentum space since it can 
be expressed as a spherical harmonic multiplied by a radial integral |131l 
Eqs. (4.3) and (4.4)]: 

F^(P) = fMY^{p/p), (4.6a) 

POO 

hip) = (~iY 2 p- 1/2 / r 3 / 2 J e2+1/2 (pr)f h (r)dr. (4.6b) 
J 

In the same way, we obtain the following Hankel-type integral representation 
for the radial part f &2 (r) of F™ 2 (r) [TED Eq. (4.5)]: 

poo 

fe 2 (r) = i^r- 1 ' 2 / p 3/2 Je 2+1/2 (rp) f\(p) dp . (4.7) 
Jo 
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Let us now consider the product y™ 1 (V)i 7 ^ 12 (r). With the help of (gUjl 
and (|A.3() . we obtain the following Fourier integral representation: 

y^(V)F^{r) = i^(2vr)- 3 / 2 J e irp y™ 1 (p) F™ 2 (p) d 3 p . (4.8) 

Next, we use the Rayleigh expansion Q4.5JI and replace the spherical harmon- 
ics by Gaunt coefficients according to (|C.4|) . This yields |1311 Eq. (4.7)]: 

yj?(y)F™{r) 

^max 

= J2 (2) ( im i + m 2 \hm 1 \l 2 m2) Y £ mi+m2 (r/r) 
f—t ■ 

/•OO 

x r- 1 ' 2 / p e ^ 2 J t+1/2 (rp) hM dp . (4.9) 



Comparison of (|2.3|) and (|4.9|) yields the following integral representation for 
the radial function 7^ 2 (r) [Ml Eq. (4.8)]: 

/■OO 

7UW = i^r-V2 / /'+ 3 / 2 J, +1/2 M/ <2 (p)dp. (4.10) 
Jo 

Now, all we need is a differential operator in r that generates the integral 
in ()4.1Uj) from the integral in (|4.7|) , With the help of known properties of 
the Bessel function J v {z) or by direct differentiation techniques (see j!31t 
Section III]), this can be accomplished relatively easily |131l Eqs. (3.29), 
(4.15) - (4.18), and (4.24)]: 

7 | i|a(r) = y {-MU-a{i)-l/2) q 2q/l+fo _ 2q 

g=0 ^ 



1 d V 1 " 5 fi 2 (r 



T d.7* / y***2 



(4.11) 



r I -— I r 

r 



dr ) 



(i 


-) 


V r 


dr / 


d 




dr 






1 d 


'( 


r dr 



^^^Mri) ^ +1 /*(r) (4.13) 

£— 1 | | ^1-^2+3^+1 
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1 f--Y"' 

\r dr J 



X^^^J ^ + V. 2 (r) (4.14) 



r dr ) \ r dr 

< r 2 ' + 1 ('^V (4.15) 
r dr / r t2 



M 2 

£ 

s=0 



X 



-A^UAh + 1/2). rrei ^ 2s ^ 



ld_ 

r dr 



*l—s 



r i2+1 fi 2 (r). (4.16) 



The abbreviations Al, Ah, Ah, and a(£) are defined in ()G6|) - (|(L9j) . Obvi- 
ously, the representations ()4.14|) and Q4.15|) make sense only if either l 2 > £ 
or I > l 2 hold. 

The expressions Q4.11JI - (|4. 16f) for 7^ 2 (r) all look quite different. Never- 
theless, their equivalence can be shown explicitly with the help of summation 
theorems for generalized hypergeometric series with unit argument (see the 
discussion in connection with Eqs. (4.19) - (4.24)]). 

Alternative expressions for jf g (r) as well as for more general radial 
functions, that occur in the case of the product V 2n y™ 2 (V)F £ ™ 2 (r) with 
n £ No, were considered by Santos [SHj, Bayman Stuart |112j . Niukkanen 
[70] . Weniger and Steinborn jl.Slj . and Rashid 

The expressions (|4.11j) - 1)4. 16j) for 7^ 2 (r) all have a manageable com- 
plexity and are well suited for practical applications. Nevertheless, it is more 
convenient to compute the product 3^™ 2 (V)i^ m2 (r) or also the more general 
product V 2n 3^ 2 (V)F™ 2 (r) with n e N via j33> if the function $ £a (r) 
defined in ()3.5|) is easily accessible and has a sufficiently simple form as in 
the case of the irregular solid harmonic and the Gaussian function according 
to ()3.9j) and (|3.11|) . respectively, or in the case of B functions according to 

{E3J). 

With the help of Fourier transformation, it is easy to obtain the explicit 
expression for the product yj l iy)f{r)g{r) which could be called the Leibniz 
theorem of the spherical tensor gradient operator and which was originally 
derived by Dunlap [22J Eq. (10)], albeit in a somewhat cryptic way. My 
derivation is inspired by Grotendorst and Steinborn |39l Appendix A] who 
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had derived a more general expression involving spherical tensor gradient 
operators in connection with the Fourier transform of a two-center density. 

For that purpose, let us express both / and g as inverse Fourier integrals 
according to (|A.3|) and use (|4.4I) : 

^ m (V)/(r)<7(r) = J J e^^f(p)g(p)d 3 pd 3 q (4.17) 

= (2^r 3 | J e ir -( p+<J )^ m (i[P + g])/»5(p)d 3 pd 3 q. (4.18) 

The dependence of the regular spherical harmonic on the two momentum 
variables p and q can be decoupled with the help of the well known addition 
theorem of the regular solid harmonics (a particularly simple derivation of 
this addition theorem based on the the spherical tensor gradient operator 
can be found in 120, Section 6]): 

£r Q (l/2jA+l(l/2j£-A+l 
A 

x <MA -ii\t- Am + ^^WO'). (4.19) 

n=—\ 

The Gaunt coefficients in this addition theorem can be expressed in closed 
form (compare for instance |120| Eq. (6.5)] which unfortunately contains a 
typographical error: The factor (2^2 — i\ + 1) in the numerator of the square 
root on the right-hand side has to be replaced by (2^2 — 2£i + 1).) In this 
way, we obtain Steinborn's factor-free version of the addition theorem of the 
regular solid harmonics (1101 Eq. (9)]. 

If we now insert (|4.19j) into 1)4. 18|) and use (|4.4|l . we obtain: 

x(2vr)- 3 / 2 | e ir Py^(ip)f(p)d 3 p 

x (2vr)- 3 / 2 J e-* ypfiiq) g(q) d 3 q (4.20) 
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aT (V2)a+i(1/2)/-a+i t^A 



/!= — A 

(2vr)- 3 / 2 ^(V) / e- p /»d 3 p 



(4.21) 



x [^"(V) /(r)] [3K"(V) fl(r)] • (4.22) 

Thus, the Leibniz theorem ()4.22|) of the spherical tensor gradient operator 
is nothing but the addition theorem (|4.19|) of the regular solid harmonic in 
momentum space. 



5 Reduced Bessel Functions 

In Section|31it was shown that the application of the spherical tensor gradient 
operator ^^"(V) to the Coulomb potential or to a scalar Gaussian yields 
remarkably compact results according to (|3.9|) or (|3.11|) . respectively. 

Other functions, to which the spherical tensor gradient operator can 
be applied with remarkable the so-called reduced Bessel functions, 

whose use in quantum chemistry had been proposed by Shavitt jlUll Eq. (55) 
on p. 15], and their anisotropic generalizations, the so-called B functions, 
which had been introduced by Filter and Steinborn 30, Eq. (2.14)]. Detailed 
discussions of the mathematical properties of reduced Bessel functions and 
B functions can be found in my PhD thesis |117j and in the PhD thesis of 
Homeier [44] . 

If K v (z) is a modified Bessel function of the second kind |61| p. 66], the 
reduced Bessel function is defined by |1U31 Eqs. (3.1) and (3.2)] 

k u (z) = (2/vr) 1 / 2 z v K v (z) . (5.1) 

If the order v of a reduced Bessel function is half-integral, v = n + 1/2 
with n £ No, the reduced Bessel function can be written as an exponential 
multiplied by a terminating confluent hyper geometric series \F\ |13L)l Eq. 
(3.7)]: 

1+1/2(2) = 2™(l/2) n e^ 1 F 1 (-n;-2n;2z). (5.2) 
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The polynomial part in (|5.2[) was also treated independently in the math- 
ematical literature, where the notation Q n (z) = e z k n +i/2( z ) is used. To- 
gether with some other, closely related polynomials, the @ n (z) are called 
Bessel polynomials [SHI- According to Grosswald [SEj) they are applied in 
such diverse fields as number theory, statistics, and the analysis of complex 
electrical networks. 

The Bessel polynomials ® n (z) occur also in a completely different math- 
ematical context: In the book by Baker and Graves- Morris p. 8] on Pade 
approximants, it is remarked that Pade had shown in his seminal thesis |81j 
that the Pade approximant [n/m] to the exponential function exp(z) can be 
expressed as the ratio of two terminating confluent hypergeometric series 
Eq. (2.12)]: 

[n/m] = -— ; r, n,mGN . (5.3) 

\t\\—n\—n — m\—z) 

Accordingly, the diagonal Pade approximant with n = rato the exponential 
function can be expressed as the ratio of two Bessel polynomials: 

As an anisotropic generalization of the reduced Bessel function with half- 
integral order, the so-called B function was introduced by Filter and Stein- 
born 30, Eq. (2.14)], 

B™ e (a,r) = [2 n+< (n + *) ! ]~ 1 fc n -i/ 2 (<*r) ^ m (ar) , a>0,n£Z. (5.5) 

B functions are a fairly large class of exponentially decaying functions. 
For n G N, they are suited to serve as trial functions in LCAO-MO calcu- 
lations. Since, however, B functions have a much more complicated math- 
ematical structure than for example Slater- type functions, whose molecular 
multicenter integrals functions are notoriously difficult, it is by no means 
obvious that anything can be gained by using B functions instead of Slater- 
type functions as basis functions. However, B functions possess a Fourier 
transform of remarkable simplicity: 

B™ e ( a ,p) = (2vr)- 3 / 2 | e^ r B^ e ( a ,r)d 3 r 

= ( 2 M 1/2 [a 2 +P 2 ]n+w yr(-ip)- (5.6) 
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This is most likely the most consequential and certainly the most often 
cited result of my PhD thesis (1171 Eq. (7.1-6) on p. 160]. Later, (|5.6|) was 
published in 129, Eq. (3.7)]. Independently and almost simultaneously, the 
Fourier transform 1)5.6)1 was also derived by Niukkanen |72| . 

The exceptionally simple Fourier transform 1)5.6)) gives B functions a 
unique position among exponentially decaying basis functions, and it also 
explains why other exponentially decaying functions like Slater-type func- 
tions, bound state hydrogen eigenfunctions or other functions sets based on 
generalized Laguerre polynomials can all be expressed in terms of finite lin- 
ear combinations of B functions (details and further references can found in 
[TT8l Section IV] or [1221 Section 4]). 

As discussed in Section 0] Fourier transformation is one of the principal 
approaches for the treatment of multicenter integrals )82) 137) . 

Thus, the simplicity of the Fourier transform ()5.6)) makes it plausible 
that multicenter integrals of B functions can normally be evaluated more 
easily than the analogous integrals of other exponentially decaying functions 
as for instance Slater-type functions, For example, by inserting ()5.6)) into 
the inverse Fourier integral on the right-hand side of 1)4.1)) . we immediately 
obtain the following, extremely simple expression for the convolution integral 
of two B functions with equal scaling parameters: 



This expression was originally derived by Filter and Steinborn |291 Eq. (4.1)] 
with the help of an addition theorem. The summation limits £ m i n and £ max 
are given in l)C.5|) . and A£ is defined by (|C6)) . 

In recent years, some significant progress has been achieved with respect 
to molecular multicenter integrals of B functions (see for example the ar- 
ticles by Steinborn, Homeier, Fernandez Rico, Ema, Lopez, and Ramirez 
)1U5) . and Steinborn, Homeier, Ema, Lopez, and Ramirez jlU4j ). Particu- 
larly promising seems to be the approach of Safouhi who - starting from the 
Fourier transform ()5.6)) - converts complicated multicenter integrals of B or 
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Slater-type functions to multi-dimensional integral representations involving 
nonphysical variables that have to be evaluated by numerical quadrature. At 
first sight, this does not look like a good idea because of the oscillatory nature 
of the multi-dimensional integral representations, which makes the straight- 
forward application of conventional quadrature methods difficult. However, 
these computational problems can be overcome if the quadrature schemes are 
combined with suitable nonlinear extrapolation methods. Based on previous 
work of SidipL02 and of Levin and Sidi [SU] on extrapolation methods for 
numerical quadrature schemes, Safouhi succeeded in developing some extrap- 
olation techniques specially suited for his needs. This permits a remarkably 
efficient and reliable evaluation of complicated molecular multicenter inte- 
grals via oscillatory (Fourier based) integral representations (see for example 
the recent articles by Berlu and Safouhi 0121 El, Berlu, Safouhi and Hog- 
gan HU, Safouhi jSHl EH EH EH E23 EI] , Safouhi and Hoggan (HI ESI EU 
and references therein). Safouhi's work can also be viewed as a convincing 
demonstration of the practical usefulness of extrapolation and convergence 
acceleration techniques in molecular electronic structure theory. 

It follows at once from (|5.6|) that a B function can be expressed as an 
inverse Fourier integral according to 

B»,r) = / dV (5.8) 



27T 2 J [ a 2 +p 2]n+e+l 

Comparison of Q4.4|) and ([5.8)1 shows that the application of the spherical 
tensor gradient operator to a scalar B function yields a nonscalar B function 
[TT71 Eq. (7.1-10) on p. 161]: 

B™/a, r) = (47^/2 (_ a )-< y™(y) B° n+i>0 (a, r) . (5.9) 

This as well as several other related results can also be deduced from known 
properties of the modified Bessel function K u (z) via ([3.4)1 (see for example 
H2SI Eq. (4.12)]): 

It is also quite easy to derive an explicit expression for the product 
3^^ ll (V)i?™ 2 £2 (r) via ([3.8|) since the application of higher powers of the 
Laplacian V 2 to a B function poses no problems. It follows at once from the 
integral representation 1)5. 8 Jl that the differential operator 1 — a~ 2 V 2 of the 
modified Helmholtz equation functions as a ladder operator [1311 Eq. (5.6)]: 

[1 - a" 2 V 2 ] B£(a, r) = B^/a, r) . (5.10) 
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As discussed in more details in Section H3 this relationship holds also if the 
indices n and £ satisfies n + 1 < 0, i.e., for B functions that are derivatives 
of the three-dimensional delta function according to (|6.19|) . 

Thus, the binomial expansion of a~ 2u V 2l/ in powers of 1 — a~ 2 V 2 in 
combination with (|5.1U|) yields |131l Eq. (5.7)]: 

a -2u V 2u K ^ r) = (-1)* ( v \ B™_ t/ (a, r) . (5.11) 

t=o ^ ' 

If we now combine (|3.8jl with ()5.11|) . we immediately obtain the following 
compact linear combination of Gaunt coefficients, binomial coefficients, and 
B functions [HD Eq. (6.25)]: 

^max 

^•(V)B; 2 ya,r) = (-a)* £ ^ {irm + m 2 \l x m x \i 2 m 2 ) 



£i ma 

x E ( t) CVtV t >,r) . (5.12) 



It follows from the summation limits (|C.5|) that A£ defined in (|C.6|) is either 
a positive integer or zero. 

As already remarked above, all the commonly used exponentially decay- 
ing functions like Slater-type functions, bound state hydrogen eigenfunctions 
or other functions sets based on generalized Laguerre polynomials can all be 
expressed in terms of finite linear combinations of B functions (see for ex- 
ample P21 Section IV] or P22 Section 4]). Thus, it follows from (t5~T2l that 
the product of a spherical tensor gradient operator and one of these expo- 
nentially decaying functions can be expressed as a finite linear combination 
of B functions. 



6 Spherical Delta Functions 

Classically, the domain of the spherical tensor gradient operator consists 
of the differ entiable functions / : M 3 — > C, although we are in practice 
only interested in differentiable irreducible spherical tensors of the type of 
1)2.2(1 . However, as for example argued in Dirac's classic book §15], the 
functions defined in the sense of classical analysis do not suffice in quantum 
theory. It is necessary to use also more general mathematical objects, the 
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so-called generalized functions or distributions, whose mathematical theory 
was rigorously formulated by Schwartz (see and references therein). 

The best known nonclassical generalization of a function / : M 3 — > C is 
the three-dimensional Dirac delta function that can be defined as a general- 
ized solution of the Poisson equation of a unit point charge |531 Eq. (1.31)]: 

V 2 - = -4n6(r). (6.1) 

r 

This Poisson equation also expresses the well known fact that the Coulomb 
potential is the Green's function of the three-dimensional Laplace equation. 

The prototype of a distribution, which is also an irreducible spherical 
tensor of rank £, is the so-called spherical delta function (see for example 
[HS1 Eq. (30)]): 

5T(r) = {2 {~X)l\ yr{V)6{r) - (6 ' 2) 

The spherical delta function can also be obtained by applying the Laplacian 
to an irregular solid harmonics (see for example |85l Eq. (29)]): 

V 2 Zr(r) = -^5T(r). (6.3) 

This follows at once from (fo~T|) . (l6~2]l . and the fact that V 2 and ^ m (V) 

commute. Comparison of (|6,lj) and (|6.3|) shows that the spherical delta 
function can be viewed as a generalized solution of the Poisson equation of 
a unit multipole charge. 

The properties of generalized functions of the type of the spherical delta 
function can be understood most easily with the help of Fourier transforma- 
tion. As already mentioned in Section |IJ Fourier transformation is defined 
in the sense of classical analysis only for functions that are absolutely inte- 
grable, i.e., which belong to the function space L 1 (M 3 ). By means of a suit- 
able limiting procedure, Fourier transformation can be extended uniquely 
to give a unitary map from the Hilbert space L 2 (M 3 ) of square integrable 
functions onto itself (see for example |841 Theorem IX. 6 on p. 10]). A further 
extension of Fourier transformation to the space of tempered distributions 
is also possible (see for example |841 Theorem IX. 2 on p. 5]). 

The extensibility of Fourier transformation to nonclassical function spaces 
is very important for our purposes. For example, the Coulomb potential is 
neither absolutely integrable nor square integrable. Nevertheless, it is pos- 
sible to define its Fourier transform in the sense of distributions |36l Eq. (2) 
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on p. 194]: 

(2/vr) 1 /2 „,„ r 



(2vr)- 3 / 2 / d 3 r. (6.4) 



p2 

This relationship makes it possible to convert multicenter integrals, which 
describe the Coulomb interaction of classical or nonclassical charge densities, 
into momentum space integrals according to (|4.3|) . This is a common ap- 
proach in the case of exponentially decaying basis functions (see for example 
126 or |4L)| and references therein). 

In the same way, we obtain the Fourier transform of the irregular solid 
harmonic, which again holds in the sense of distributions [1261 Eq. (3.11)]: 



Z?(p) = (2vr)- 3 / 2 J e- ip r Z^(r)d 3 i 



It is also quite instructive to study the Fourier transforms of distributions 
as for instance the three-dimensional delta function or the spherical delta 
function. If we set f(r — r') = 5{r — r') in the convolution integral l|4.1jl . we 
see that the Fourier transform of the delta function is a constant: 

5{p) = (2tt)- 3 / 2 J e- ip - r 5{r)d 3 r = (2tt)- 3/2 . (6.6) 

Similarly, we find that the Fourier transform of the spherical delta function 
is essentially a regular solid harmonic in momentum space: 

In physics, it is common to introduce the three-dimensional delta function 
via the Poisson equation (|6.1|) . Both the Poisson equation (|6,lj) as well as 
its anisotropic generalization l)6.3j) follow from the Fourier transforms (|6.4j) 
and (|6.5|) . respectively. We only have to take into account that the function 
l/p 2 , which occurs in both (|6.4j) and (|6.5[) . is canceled by the Laplacian V 2 , 
whose Fourier transform is —p 2 . 

It is, however, just as well possible to introduce the there-dimensional 
delta function as well as the spherical delta function via the differential oper- 
ator 1 — a _2 V 2 of the modified Helmholtz equation, whose Fourier transform 
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is given by [a 2 +p 2 ]/a 2 . This follows at once from that fact that the Yukawa 
potential |14Uj . which is an exponentially screened Coulomb potential with 
screening parameter a, is also a special B function according to 

= ak_ 1/2 (ar) = (4vr) 1/2 a B$ fi (a, r) . (6.8) 



Obviously, the Yukawa potential approaches the Coulomb potential in the 
limit of vanishing screening, i.e., for a — ► 0. 

If we set n = £ = m = 0m (J5.6|) . we essentially obtain the Fourier 
transform of the Yukawa potential: 

B%(a,p) = (27T)- 3 / 2 | e- i ^ J B ° (a,r)d 3 r 

= (2vr 2 )-^ . (6.9) 
[or +p \ 

Next, we use (|4.1j) to compute the convolution of a relatively arbitrary func- 
tion / : R 3 -> C with the Yukawa potential Q3D Eq. (6.9)]: 

J B%(a,[r-r'])f(r')dr = J B° ofi (a,p) f(p) d 3 p 

= w^IJr? f ~^ p - (6 - 10) 

If we now apply the differential operator 1 — a _2 V 2 to the convolution inte- 
gral, interchange integration and differentiation, and use ()5.1Uj) . we see that 

(a, r) is proportional to the three-dimensional 
delta function Q3D Eq. (6.10)]: 

[1-«" 2 V 2 ] f B%(a,[r-r f ])f(r')dr 

= J B°_ 1;0 (a,[r-r'])f(r')dr 

= ^r 3/2 1 * irp m i a p = ^ /« • 

cr J or 

We thus obtain the following exponentially screened variant of the Poisson 
equation (|6.1|) : 

[1 - a" 2 V 2 ] B Q Q< (a, r) = <S(r) . (6.12) 
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This relationship can also be expressed in terms of the Yukawa potential: 

— OiT 

[a 2 -V 2 ]^— = 4vr5(r). (6.13) 

For o = 0, we obtain the Poisson equation ()6.1|) . 

The screened Poisson equation (|6,13|) expresses the well known fact that 
the Yukawa potential is the Green's function of the modified Helmholtz equa- 
tion (see for example 2, Table 16.1 on p. 912]). 

An anisotropic generalization of the approach described above is also 
possible. If we set n = — I in 1)5. 9|) . we see that the application of the 
spherical tensor gradient operator to the Yukawa potential yields the so- 
called modified Helmholtz harmonic |133| Eq. (6.9)]: 

B m e/ (», r) = (4n)^ 2 (- a y e ^(V) B° ofi (a, r) (6.14) 

In the limit of vanishing screening, the modified Helmholtz harmonic ap- 
proaches an irregular solid harmonic according to |126| Eq. (3.10)] 

Zf{v) = [{21 - l)!!]- 1 lim [a^B m ^a,r)] . (6.15) 

If we set n = — t in ()5.6|) . we obtain the Fourier transform of the modified 
Helmholtz harmonic |129| Eq. (A.l)]: 

B^{a,p) = {2/v) l l \ a ~*~ l y?{-ip). (6.16) 

[or + p z \ 

Both the Fourier integral producing this relationship as well as the inverse 
Fourier integral representation for B™p Ja, r) do not exist in the sense of clas- 
sical analysis. In [1291 Appendix] it was, however, shown that the classically 
divergent Fourier integral representation for B™» Ja,r) can be regularized 
by employing a suitable rational cutoff function. 

If we now proceed as in the case of (|6.11jl and also use (|4.4|) . we see that 
B r ^ e __ l e (a, r) = [1 — a~ 2 V 2 ]B m £ e (a, r) is proportional to the spherical delta 
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function HSU Eq. (6.17) - (6.19)]: 

[l_ a -2 V 2] J Bi e/ (a,[r-r'])f(r>)dr 



Bi e _ 1>e (a,[r-r'})f(r')dr 
(2Ar)V2 



a e + 3 



J jr- PyT{ _ ip) f - {p)d 3 p 



47T 



(-ir^sOT/M. (6.i7) 



c^+ 3 

We thus obtain the following exponentially screened variant of the anisotropic 
Poisson equation (|6,3j) : 

[l-a" 2 V 2 ] B™ e/ (a,r) 

Air Att 

= (-l)'^T3^ m (V)^(r) = (2e-l)\\5T(r). (6.18) 

In view of this relationship and also because of 1)5. 1U() it makes sense to 
define a distributional B function as the following derivative of the three- 
dimensional Dirac delta function |131| Eq. (6.20)]: 

B m k _ E/ (a, r) = (2£ ~ 1 + ) 3 !47r [1 - a" 2 V 2 ] ^ «5f (r) , kN. (6.19) 

The fact that the distributional B function -B^_ 1 » is proportional to the 
spherical delta function can also be seen by setting ri\ = —l\ — 1 in the 
convolution integral (|5.7|) . Then, we obtain apart from a different prefactor 
the expression (I5TT2]) for the product ^(V) B™ 2 4 (a,r). 

The distributional nature of some B functions follows also from the fact 
that the Fourier transforms (|5.6[) of B functions satisfy for all n £ Z and for 
all n € No the following functional equations |1311 Eqs. (6.21) - (6.23)]: 

B™ E (a,p) = [a 2 /(a 2 +p 2 )]B^ e (a,p), (6.20a) 
B™ e (a,p) = [(4n)^/a e ]yr(-ip)B° n+efi (a,p), (6.20b) 
B°_ lfi (a,p) = a' 3 (2ir 2 )- 1/2 . (6.20c) 

These functional equations in momentum space show that there is an in- 
timate relationship between B functions, the differential operator of the 



Ernst Joachim Weniger: The Spherical Tensor Gradient Operator 



7: Addition Theorems 



29 



modified Helmholtz equation, the spherical tensor gradient operator, and 
the three-dimensional delta function. 

There is a simple and intuitive interpretation of distributional B func- 
tions. Because of the factorial (n + £)\ in the denominator on the right-hand 
side of (|5.5[) . B functions are defined in the sense of classical analysis only if 
n + £ > holds. However, the definition of a B function remains meaningful 
even for n + £ < 0. If r ^ 0, the value of B™ k _ i£ (r) with k € N is because 
of the singular factorial (—k)l = T(—k + 1) zero, but for r = 0, its value 
is oo/oo and therefore undefined. Thus, the mathematical nature of a B 
function with n + 1 < as well as its value for r = has to be analyzed with 
the help of suitable limiting procedures. 

7 Addition Theorems 

In many subfields of physics and physical chemistry - for example in electro- 
dynamics _53 , in classical field theory , or in the theory of intermolecular 
forces jl 11] - an essential step towards a solution of the problem under con- 
sideration consists in a separation of variables. 

Principal tools, which can accomplish such a separation of variables, 
are so-called addition theorems. These are expansions of a given function 
f(r ± r') with r, r' € M 3 in products of other functions that only depend on 
either r or r', 

In view of the importance of addition theorems, it is not surprising that 
there is a very extensive literature. Consequently, any attempt to provide 
a reasonably complete bibliography would clearly be beyond the scope of 
this article. The interested reader is referred to the long, but nevertheless 
incomplete lists of references in |12U1 I122| . 

Addition theorems have played a major role in my own research. I first 
applied addition theorems for the evaluation of some multicenter integrals of 
B functions in my Master's thesis |116j . which was published in condensed 
form in |l()7j . In my PhD thesis |117j and also afterwards, I preferred Fourier 
transformation for the evaluation of multicenter integrals of B functions, but 
I later worked on the derivation of addition theorems |4fi[ I118| I12()| 11221 M'A'M 

H5B| . 

In atomic or molecular calculations, we are predominantly interested in 
irreducible spherical tensors of the type of (|2.2|) . Moreover, the convenient 
orthonormality of the spherical harmonics makes it highly desirable that the 
functions of either r or r', which occur in the expansion of f(r ±r'), are also 
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irreducible spherical tensors. Thus, the addition theorems, we are interested 
in, are expansions in terms of spherical harmonics with arguments 8,4> = r/r 
and 6',4>' = r'/r', respectively. 

The best known example of such an addition theorem is the Laplace 
expansion of the Coulomb potential in terms of spherical harmonics: 

1 1 A=0 fi=-X 

r < = min(r, r ) , r> = max(r, r') . (7.1) 

The Laplace expansion leads to a separation of the variables r and r' . How- 
ever, its right-hand side depends on r and r' only indirectly via the vectors 
r< and r> which satisfy |r<| < |r>|. Hence, the Laplace expansion has a 
two-range form, depending on the relative size of r and r'. This is a com- 
plication which occurs frequently among addition theorems. As discussed 
in more details in |12Ul I122| . addition theorems have a two-range form if 
they are pointwise convergent three-dimensional Taylor expansions and if 
the function f(r±r'), which is to be expanded, is not analytic at the origin. 

The undeniably troublesome two-range form of an addition theorem can 
be avoided if / : M 3 — » C belongs to the Hilbert space L 2 (R 3 ) of square in- 
tegrable functions or to other, closely related function spaces as for example 
Sobolev spaces that are proper subspaces of L 2 (IR 3 ) (compare for instance 
|31| 1401 1118[ I132| and references therein) . For the sake of simplicity, let us 
assume that a discrete function set {^^(f )}n,4m is complete and orthonor- 
mal in the Hilbert space L 2 (M 3 ). Then, an addition theorem for f(r ± r'), 
which converges in the mean with respect to the norm of L 2 (M 3 ), can be 
derived by expanding / in terms of the functions £(?*)}n,£,m : 

f(r ± O = E C ™Af-> r ') *W . ( 7 - 2a ) 

nlm 

C^(/;r') = J [¥&(!■)]• f(r±r')d 3 r. (7.2b) 

Such an expansion is a one-range addition theorem since the dependence on 
r is entirely contained in the functions ^^Ar), whereas r' occurs only in the 
expansion coefficients C™Jf; r') which are overlap integrals. With minimal 
modifications, this approach works also if the functions ^(»*)}n,^m are 
complete and orthonormal in a suitable Sobolev space. 
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At first sight, it looks as if one-range addition theorems of the type of 
(|7.2|) are clearly superior to two-range addition theorems, but a balanced 
assessment of their relative merits is not so easy. Firstly, one-range addi- 
tion theorems usually have a more complicated structure than two-range 
addition theorems (typically, one-range addition theorems contain one addi- 
tional infinite summation). Secondly, one-range addition theorems normally 
converge only in the mean with respect to the norm of the Hilbert space 
L 2 (M 3 ), but not necessarily pointwise. In some cases, this can lead to con- 
vergence problems. The probably the most severe disadvantage is that the 
approach sketched above works only if / is an element of a suitable Hilbert 
or Sobolev space. Many functions of considerable practical importance such 
as the Coulomb potential or the irregular solid harmonic do not belong to 
L 2 (M 3 ), let alone to a suitable Sobolev space. Therefore, it is not possible to 
derive one-range addition theorems by expanding them in terms of function 
sets that are complete and orthonormal with respect to a scalar product that 
involves an integration over the whole M 3 . 

In this article, I will only discuss addition theorems that converge point- 
wise, i.e., which can be viewed to be three-dimensional Taylor expansions. 
Moreover, I will concentrate on addition theorems that can be derived with 
the help of the spherical tensor gradient operator y7 l (V). 

As discussed in Section El the irregular solid harmonic Z™(r) can ac- 
cording to ()3.9|) be obtained by applying y™(V) to the Coulomb potential. 
Thus, it should be possible to derive the addition theorem of the irregu- 
lar solid harmonic by applying either 3^™(V<) or y™(S7>) to the Laplace 
expansion 1)7.1)) . As shown in |133) Section IV], this is indeed possible. More- 
over, it was shown in (1331 Sections V and VI] that the addition theorems 
of the Helmholtz and the modified Helmholtz harmonics can also be derived 
by applying either 3^™(V '<) or y™(y>) to the simpler addition theorems of 
the corresponding scalar functions. It is also possible to derive the addition 
theorem of B functions by applying the spherical tensor gradient operator 
to the addition theorem of the reduced Bessel functions [TJjBJ Section V] . 

The approach described above requires that for a given irreducible tensor 
a scalar function can be found which satisfies ((3.5)) and which possesses a 
suitable addition theorem. This obviously limits the practical usefulness 
of this approach. It is, however, possible to derive addition theorems of 
essentially arbitrary irreducible spherical tensors from the scratch with the 
help of the spherical tensor gradient operator. 

As is well known, addition theorems can formally be obtained by doing 
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three-dimensional Taylor expansions (see for example [I'll p. 181]): 

f(r + r') = Yl ^—^Kr) = ^/(r). (7.3) 

n=0 

Thus, the translation operator 

g r' V _ e x'd/dx e y'd/dy e z'd/dz {7 A) 

generates f(r + r') by doing a three-dimensional Taylor expansion of / 
around r with shift vector r' . Moreover, the variables r and r' are sep- 
arated. Thus, the series expansion ()7.3|) is indeed an addition theorem. 

We could also expand / around r' and use r as the shift vector. This 
would produce an addition theorem for f(r + r') in which the roles of r and 
r' are interchanged. Both approaches are mathematically legitimate and 
equivalent if / is analytic at r, r', and r + r' for essentially arbitrary vectors 
r,r' £ R 3 . Unfortunately, this is normally not true. Most functions, that 
are of interest in the context of atomic and molecular quantum mechanics, 
are either singular at the origin or are not analytic at the origin. Obvious 
examples are the Coulomb potential, which is singular at the origin, or the Is 
hydrogen eigenfunction, which possesses a cusp at the origin. In fact, all the 
commonly used exponentially decaying functions as for example Slater-type 
functions or also B functions are not analytic at the origin. 

The reason for the non-analyticity is that the three-dimensional distance 

1 /2 

r = [x 2 + y 2 + z 2 ~\ is not analytic with respect to x, y, and z at the 
origin r = 0. This implies that all odd powers r 2n+1 with n £ No are also 
not analytic at the origin (compare also the discussion related to |118| Eq. 
(5.9)]). In contrast, r 2 = x 2 + y 2 + z 2 and the regular solid harmonic y™(r) 
are analytic since they are polynomials in x, y, and z. Consequently, a Is 
Gaussian function exp(— ar 2 ), which possesses an expansion in even powers 
r 2n , is analytic at r = 0, but a Is Slater-type function exp(— ar) is not. 

Thus, for the derivation of addition theorems for functions, that are not 
analytic at the origin, we have to use the translation operator in the following 
form, 

f(r< +r>) = V lr< ' j >} /(r>) = e^ v > /(r>) , (7.5) 

n=0 

where |r<| < |r>[. In this way, the convergence of the three-dimensional 
Taylor expansion is guaranteed provided that / is analytic everywhere except 
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possibly at the origin. Thus, the non-analyticity of B functions and of all the 
other commonly occurring exponentially decaying functions at the origin has 
a far-reaching consequence: Their pointwise convergent addition theorems 
must have a two-range form in order to guarantee convergence (see also |122| 
and references therein). 

Prom a practical point of view, the translation operator e r< v> in its 
Cartesian form does not seem to be a particularly useful analytical tool. We 
are usually interested in addition theorems of irreducible spherical tensors, 
which are defined in terms of the spherical polar coordinates r, 8, and 4>. 
Differentiating them with respect to the Cartesian components x, y, and z 
would lead to extremely messy expressions and to difficult technical prob- 
lems. Thus, it is tempting, but nevertheless superficial to conclude that the 
translation operator e r< v> provides only a formal solution to the problem 
of separating the variables r and r' of a function f(r + r 1 ). 

The crucial step, which ultimately makes the Taylor expansion method 
practically useful, is the expansion of the translation operator e r< v> in 
terms of differential operators that are irreducible spherical tensors. The 
starting point is an expansion of exp(a-b) with o,b6t 3 in terms of modified 
Bessel functions and Legendre polynomials |HB p. 108]: 



e 



,a-6 ^ab cos 9 / 



1/2 00 

(■fb) E ( 2£ + !) ^+i/2(a*0 ^(cos 6) . (7.6) 
a e=o 



Next, the series expansion for the modified Bessel function Ig + i/2 EH P- 
66] is inserted into Q7.6JI . and spherical harmonics are introduced with the 
help of the so-called addition theorem of the Legendre polynomials (see for 
example j!2l p. 303]). This yields the following expansion of e a b in terms of 
regular solid harmonics and even powers of the vectors a and b: 



00 « 00 2k u2k 

e=o m=-e k=o ' y 1 ' t+K + 1 

The powers a 2k and b 2k are irreducible spherical tensors of rank zero, and 
the solid harmonics are tensors of rank I. 

The expansion (|7.7|) is obtained from (|7.6|) by rearranging the Cartesian 
components of the vectors a and b. Accordingly, it holds for essentially 
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arbitrary vectors a,b E R 3 , and we can choose a = r < and b = V>: 

_ oo £ oo ok Tj2k 

e-v> = 2»£ £ pr (, <)r 3T(V > ) £ 2<+a ,r ( l/2) <+t+1 • < 7 - 8) 



1=0 m=-t k=0 



It seems that this expansion was first published by Santos [23 Eq. (A. 6)], 
who emphasized that this expansion should be useful for the derivation of 
addition theorems, but apparently never used it for that purpose. 

It follows from the expression (|B.2|) of the Laplacian V 2 in spherical polar 
coordinates and from the tensorial nature of the spherical tensor gradient 
operator (compare Q2.3|) and the numerous expressions for j% i 2 ( r ) given in 
Section^} that in (|7.8|) we only have to do differentiations with respect to the 
radial variable r>. This greatly simplifies practical applications. In |12L)j . it 
was shown that starting from the expansion (|7.8[) the addition theorems of 
the Coulomb potential, the regular and irregular solid harmonics, and the 
Yukawa potential can be derived quite easily, and in 122 , several different 
forms of the addition theorem of the B functions was derived in this way. 

As a further demonstration of the usefulness of the expansion l|7.8j) of 
the translation operator in terms of irreducible spherical tensors I will now 
derive an addition theorem of the function r u y™(r) with u € R. 

Our starting point the the following relationship which follows at once 
from (JsQJ): 

r " y ^ r) = wrhm l yr{V,r " 2 '- <7 ' 9) 

If we combine this with ()7.8j) and linearize the spherical tensor gradient 
operators according to 1)3. 7J) . we obtain: 



e 



|r<+r>|-3r(r< + r>) = ^TUWe ^ E 



t 1= mi=-h 

e 2 oo ^k 

l2 =tfn k=0 l 

x v^ +2A/a y™ +mi (V>) r v + 21 . (7.10) 

The abbreviations A£, A£\, and A^2, which will be used in the following 
formulas and which are in all cases either zero or a positive integer, are 
defined by El 
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In the next step, we use (|7.8|) once more to obtain 
y™ +mi {V)r v+2i = (-2) i2 {-£ - v/2) h r u+2£ - 2e2 y™ +mi {r) . (7.11) 
Inserting this into (|7.1()j) yields: 

ir< + -> r mr< + ->> = ni+vl2)i E jc_ K" M " 

/max 

x ^ ( 2 ) (£ 2 m + mi\£imi\£m) 

^ fc!(l/2), 1+fc+1 

x V 2fc+2A£ 2 r ^+2«-2/ a 3;™+™! (r>) (? _ 12) 



X 



Now, the only thing, that remains to be done, is the application of the 
powers V^ fc+2A£2 of the Laplacian to r^ 21 ' 2 ^ y™ +mi (r>). For that pur- 
pose, we use the following expression which can be deduced easily from 
(TP) : 

d 2 21 + 2 9 



V^(r)3T(r) = % m (r) 



+ 



ij>(r) . (7.13) 



Here, ^(r) is a scalar function. If we set ip(r) = r a with a £ R, we find: 

VV^ m W = 4(-a/2)(-[ ( j + 2^ + l]/2)r' T - 2 ^ n (r). (7.14) 
Iterating this relationship n times yields: 

S/ 2n r°yp(r) = A n {-a/2) n {-[a + 2l+l\/2) n r a - 2n yr(r). (7.15) 
Thus, we obtain for the remaining differentiations: 

y2k+2M 2 r v+2t-2l 2 ym+ mi ^ r ^ _ ^k+A£ 2 

u + 2£-2£ 2 \ ( z^ + 2^ + l\ (2M-v s 



~ /A< 2 \ ^ 7 M 2 V 2 

2Al 1 + ^+l ^ ^+2A/ 1 -2*+l2^fm 1(r>) _ (?16) 
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With the help of some essentially straightforward algebra, it can be shown 
that the k summation in (|7.12|) can be expression by a Gaussian hyperge- 
ometric series 2-^1 (for its definition, see for instance |61l p. 37]), and we 
finally obtain the following addition theorem: 

a 00 h 

\r<+r > Yy?(r < + r > ) = —?—y: £K><)]* 

V ' n £ 1= mi=-<i 

x (2) (- 1 )' 2 (hm + milimnlim) 

e 2 =e™ in 

(-£-u/2) h /u-2A£ + 2\ (v-2A£ + 3 
X (3/2) £l V 2 ) Ai2 V 2 /Afa 

(IM-v -2A£ x -u-l 2^ + 3 r| 

x 2^1 i v ^— . — 2 — ; —t" ; ^ 

x r ^+ 2A <i+ 1 Z™ +mi (r>) . (7.17) 



The addition theorem (|7.17|) contains many simpler addition theorems 
as special cases. For example, if we set in (j7.17|) £ = m = 0, we obtain the 
addition theorem of the corresponding scalar function: 



OO 



+ r>|" = ATrr^Yj-l) 
x 2^1 



-"/2) A 



to" *' ^ 
2\-v -v-l 2A + 3 r\ 
~^~ , ^~' , ^2~' , ^l 

A 

x £ [W^w- ( 7 - 18 ) 

/Lt=— A 

Of course, we could also go the other way round: We could proceed as in 
|133j and construct the addition theorem (|7. 17|) by applying either 3^™(^<) 
or J7(V>) to (EH. 

If we set in (ITTTKl) v = -1, we only need (l/2) A /(3/2) A = 1/(2A + 1) as 
well as 2-^1 (A + 1/2,0; A + 3/2; r</r>) = 1 to obtain the Laplace expansion 
(|7.1j) of the Coulomb potential. 

If we set in (|7.18j) v = 2n with n G No, the A summation terminates 
after A = n because of the Pochhammer symbol (— n)\. In addition, the 
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hypergeometric series 2F1 i n 1)7.18)) terminates to become a polynomial of 
degree n — A in 7</r>. Moreover, r 2n is a polynomial in x, y, z and thus 
analytic. Accordingly, a distinction between r< and r> is not necessary and 
the addition theorem has a one-range form. 

A detailed analysis of all special cases of the addition theorem (|7,17[) 
would clearly be beyond the scope of this article (see also (1141 pp. 168- 
169]). Here, one must not forget that the emphasis of this Section is not on 
the derivation of the addition theorem (|7.17|) . Rather, it is my hope that the 
fairly effortless derivation of this addition theorem convinces the reader that 
the expansion (|7.8|) of the translation operator in terms of irreducible spher- 
ical tensors is indeed a highly useful analytical tool. Of course, a sceptical 
reader may well argue that it is unlikely that convenient explicit expres- 
sion for the necessary radial differentiations of arbitrary order can always be 
found, as it was the case with ()7.17)) . This certainly true. But even if we 
can only do the differentiations explicitly up to a finite maximum order, the 
tensorial expansion (|7.8|) permits at least the construction of approximations 
to addition theorems. Computer algebra systems like Maple or Mathematica 
should be helpful in this respect. 

8 Summary and Conclusions 

The regular solid harmonic y™{r) = r e Y^ n (9, 4>) is according to (|B,12|) a ho- 
mogeneous polynomial of degree i in the Cartesian components of r. Thus, 
it makes sense to define the differential operator y™(V) via (|2.1j) . i.e., by 
replacing the Cartesian components of r in ()B.12|) by the Cartesian compo- 
nents of V. 

The spherical tensor gradient operator 3^""(V) is an irreducible spherical 
tensor of rank £. This is a very consequential fact. Firstly, its application 
to a scalar function, which is an irreducible spherical tensor of rank 0, must 
produce an irreducible spherical tensor of rank £. This follows at once from 
the simplified version 1)3.4)1 of Hobson's theorem which is discussed in Section 
121 With the help of Hobson's theorem, it is also possible to derive the com- 
pact explicit expression for the product y™ 1 (V)-P^ 2 (r), where F^ 2 
is an irreducible spherical tensor with nonzero rank £2 of the type of (|2.2j) 
that also satisfies 1)3.4)1 . Secondly, the structure of products of the type of 
y™ 1 (V)F^ 2 (r) according to 1)2.3)1 can be understood completely on the basis 
of the usual coupling rules of angular momentum theory. Accordingly, the 
angular parts of such a product can be expressed in terms of Gaunt coeffi- 
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cients and spherical harmonics. Moreover, the radial parts of y™ 1 (V)F^ 2 (r) 
and of Fj£ 2 (r) are connected by relationships which only involve differenti- 
ations with respect to the radial variable r. It is this property which makes 
the spherical tensor gradient operator a practically useful analytical tool. 

Fourier transformation is one of the principal techniques for the evalu- 
ation of molecular multicenter integrals. It is also extremely useful in con- 
nection with the spherical tensor gradient operator. Under Fourier transfor- 
mation, the differential operator y™ 2 (V) produces a regular solid harmonic 
y™ 2 (ip) in momentum space. In this way, we can easily understand the 
tensorial nature of y™ 2 (S/). Moreover, many analytical manipulations can 
be done more conveniently in the momentum than in the coordinate repre- 
sentation. 

There are some scalar functions of physical interest which produce par- 
ticularly compact results if y™ 2 (V) is applied to them. Examples are the 
Coulomb potential, which produces an irregular solid harmonic according to 
()3.9)) . or the Is Gaussian function, which produces a spherical Gaussian func- 
tion according to 1)3.11)1 . Another class of scalar functions, to which y^iV) 
can be applied with remarkable ease, are the so-called reduced Bessel func- 
tions defined in (|5.1f) . and their anisotropic generalizations, the so-called B 
functions defined in 1)5.5)1 . These functions, which have played a major role 
in my own research, are discussed in Section |SJ Application of the spherical 
tensor gradient operator to a scalar B function simply produces a nonscalar 
B function according to (|5.9|) . and the product y™ 1 (V)B™ 2 i2 (a, r) can ac- 
cording to 1)5.12)1 be expressed by a simple linear combination of B functions. 
The simplicity of both ()5.9)1 and (|5.12|l follows directly from the very sim- 
ple Fourier transform 1)5.6)1 of B functions, which gives them an exceptional 
position among exponentially decaying basis functions. 

Classically, the domain of the spherical tensor gradient operator consists 
of the differentiable functions / : M 3 — > C. However, as discussed in Section 
El it makes sense to apply y™(V) also to distribution or generalized func- 
tions. This produces mathematical objects like the spherical delta function 
<5™(r), which is defined by 1)6.2)1 and which can be viewed as a generalized 
solution of the Poisson equation 1)6.3)1 of a unit multipole charge. With the 
help of Fourier transformation, the nature of the spherical delta function 
can be made transparent, and we can also understand easily why the Pois- 
son equation 1)6.1)1 holds, or - to put it differently - why convolution with 
the Coulomb potential and application of the Laplace operator are inverse 
operations. An essentially identical approach works also in the case of B 
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functions. It follows from (|5.1(Jj) or also from the functional equations (|6.2Uj) 
that the differential operator 1 — a~ 2 V 2 of the modified Helmholtz equation 
functions can be viewed as a kind of lowering operator for B functions. Since 
the Yukawa potential exp(— ar)/r is according to lj6.8|) proportional to the B 
function B® (a, r), it turns out that convolution with the Yukawa potential 
and the application of the differential operator 1 — a _2 V 2 are according to 
()6.11j) inverse operations. 

In many subfields of physics and physical chemistry, an essential step 
towards a solution of the problem under consideration consists in a sepa- 
ration of variables. Addition theorems, which are discussed in Section [7J 
are principal tools to accomplish such a separation of variables. As is well 
known, addition theorems can be obtained according to ()7.3j) by applying 
the translation operator e r v to a function f(r). However, in atomic and 
molecular electronic structure calculations, we are predominantly interested 
in irreducible spherical tensors of the type of (|2.2|) . and the convenient or- 
thonormality of the spherical harmonics makes it highly desirable that the 
functions, which occur in the addition theorem, are also irreducible spheri- 
cal tensors of the type of (|2.2|) . Accordingly, the translation operator in its 
Cartesian form ()7.4|) is not suited for our needs, since its use would lead to 
enormous technical problems. It is a much better idea to express the trans- 
lation operator in terms of irreducible spherical tensors according to (|7.8|) . 
The feasibility of this approach is demonstrated by deriving the addition 
theorem of the function r u y™(r) with v £ M. 

Let me conclude this article by some personal remarks. There can be no 
doubt that y^i^J) is a useful analytical tool since its application to a scalar 
function produces an irreducible spherical tensor of rank t. This is a highly 
advantageous feature, because now it is in principle sufficient to consider 
only multicenter integrals of scalar functions. Higher angular momentum 
states can be generated by differentiation with respect to the nuclear coordi- 
nates. Often, this is simpler that the direct derivation of explicit expressions 
for integrals of nonscalar functions. This approach has the additional ad- 
vantage that it facilitates the use of computer algebra systems like Maple 
or Mathematica, whose systematic utilization I wholeheartedly recommend. 
The completely symbolic treatment of complicated multicenter integrals is 
still too difficult for computer algebra systems, but they are already now 
very well suited for doing complicated differentiations. 

Nevertheless, one should never forget that a nice explicit expression for a 
multicenter integral does not necessarily permit its efficient and reliable eval- 
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uation. Already during the process of deriving an explicit expression, one 
should always take numerical aspects into account, because they ultimately 
decide whether a given formula is practically useful or not. Unfortunately, 
this is more easily said than done. Thus, the ability of skillfully manipu- 
lating complicated mathematical expressions does not guarantee success in 
the integral business. It is also necessary to have a detailed knowledge of 
numerical analysis. 

For example, during the work for my PhD thesis |117j , series expansions 
for multicenter integrals played a major role. Unfortunately, it often hap- 
pened that series expansions converged quite slowly (see for example |13U1 
Table II]). In principle, it is a fairly obvious idea to try to speed up the 
convergence of these expansions with the help of series transformations, but 
during my PhD thesis I only knew the linear series transformations described 
in the classic book by Knopp j^JJ, which did not accomplish much. At that 
time, I was completely ignorant about the more modern and more effective 
computational tools such as Pade approximants or other nonlinear transfor- 
mations, which often accomplish spectacular improvements of convergence. 

The situation changed radically when I did postdoctoral work with Pro- 
fessor Jiff Cfzek at the Department of Applied Mathematics of the Univer- 
sity of Waterloo. There, I worked on distributive expansions of a plane wave 
jllSj . which converge weakly with respect to the norm of the Hilbert space 
L 2 (M 3 ) or the Sobolev space W^QR 3 ). My second research topic in Water- 
loo - the summation of factorially divergent power series as they for instance 
occur as asymptotic expansions for special functions or in perturbation ex- 
pansions of quantum physics - was in the long run far more consequential, 
although I did not accomplish anything worth publishing during my stay in 
1983. 

After my return to Regensburg, I tried to apply nonlinear transfor- 
mations also to slowly convergent series expansions for multicenter inte- 
grals. In some cases, remarkable improvements of convergence were observed 

In order to understand better the power as well as the limitations of non- 
linear sequence transformations, I also worked on their theoretical properties. 
As a by-product, I was able to derive several new transformations. The ma- 
jority of these transformations was published in my long article |119| . where 
also efficient algorithms for the computation of sequence transformations are 
discussed as well as theoretical error estimates and convergence properties. 
More recent activities are described in the articles [71 11211 H231 11241 1127j . 
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Thus, it is probably justified to claim that numerical mathematics ulti- 
mately profited via cross fertilizations from the convergence problems which 
I encountered. My personal experience highlights the central importance of a 
functioning communication between mathematicians and scientists. Gener- 
alists like Professor Josef Paldus, who is at home both in mathematics and in 
physical and theoretical chemistry, help to overcome the current communica- 
tion problems and thus do an invaluable service to the scientific community. 



A Appendix: Terminology and Definitions 

For the set of positive and negative integers, I write Z = {0, ±1, ±2, . . .}, for 
the set of positive integers, I write N = {1,2,3,...}, and for the set of non- 
negative integers, I write No = {0, 1, 2, . . .}. The real and complex numbers 
and the set of three-dimensional vectors with real components are denoted 
by R, C, and M 3 , respectively. 

For the commonly occurring special functions of mathematical physics I 
use the notation of Magnus, Oberhettinger, and Soni unless explicitly 
stated otherwise. 

Double factorials are defined according to 



(2n)! 
(2n - 1)! 
0! 



2-4---2n, raeN, (A.la) 
1- 3-- -(271-1), ra€N, (A.lb) 
1!! = (-1)!! = l. (A.lc) 



Fourier transformation is used used in its symmetrical form, i.e., a func- 
tion / : M 3 — s- C and its Fourier transform / are connected by the integrals 

f(p) = (2tt)-^ 2 J e-^/(r)d 3 r, (A.2) 
f(r) = (2vr)- 3 / 2 J e ir Pf(p)d 3 p, (A.3) 

B Appendix: Spherical Harmonics 

Normally, the spherical harmonics Y™{0,4>) are introduced as the simulta- 
neous normalized eigenfunctions of the orbital angular momentum operators 

L = -M0)de sm{e) de- sTn^)^' 
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which is essentially the angular part of the three-dimensional Laplacian in 
spherical polar coordinates, 



1 d 



8 



V z = —— r z - — 

V 9 ci ' ci 9 ) 



Or 



and 



Or 




(B.2) 



(B.3) 



This determines Y^ ri (9, <fi) up to an arbitrary phase. If we for instance choose 
the phase convention of Condon and Shortley 50] Chapter III. 4], we obtain 
the following explicit expression p. 69]: 



m+\m\ 



(2l + \)(i-\m\)\ 
4tt(£+ \m\)l 



nl/2 



P l e ml (cos9)e im * . (B.4) 



Here, pJ m '(cos6>) is an associated Legendre polynomial p. 155]: 



P™{x) = (l-x 2 ) m/2 



dx 



l+m 



1 



r 2\m/2 



dx r 



Pt(x). (B.5) 



Alternative phase conventions for the spherical harmonics are discussed in 
QUI pp. 17 - 22]. 

The spherical harmonics Y™(9, (f>) are often called surface harmonics be- 
cause the angles 9 and cfi characterize a point r/r on the surface of the 
three-dimensional unit sphere. In the literature, it is common to introduce 
the so-called regular and irregular solid harmonics 



yr(r) 



- l Yn 



(B.6) 
(B.7) 



Both y™ and Z™ are defined for arbitrary vectors r € R 3 . Moreover, it is 
easy to show that the regular solid harmonics are for all r£R 3 solutions of 
the homogeneous three-dimensional Laplace equation 



v 2 /W 



o 2 



o 2 



dx 2 dy 2 dz 2 



f(r) = 0, 



(B.8) 



whereas the irregular solid harmonics are solutions for all r G R 3 \ {0}. 

In connection with the differential operator y™(V), it is more natural to 
introduce the regular solid harmonics as polynomials satisfying the Laplace 
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equation (|B.8|) . Suitable subsets of the polynomials in x, y, and z can be 
characterized by their transformation properties. For example, the homoge- 
neous polynomials Pe(x, y, z) of degree I satisfy 

Pl(rjx,r]y,r]z) = rf P e (x, y, z) , i)£C, i e N . (B.9) 

A special subset of the homogeneous polynomials of degree £ are the so-called 
harmonic polynomials Hg(x,y,z) of degree I that are also solutions of the 
homogeneous Laplace equation ()B.8|) . i.e., they satisfy 

V 2 H e {x,y,z) = 0. (B.10) 

For a given £ 6 No there are exactly 1£ + 1 linearly independent harmonic 
polynomials H^(x,y,z) (see for example 031 P- 123] or [7IJ Appendix H.3]). 
Accordingly, the regular solid harmonics span the space of harmonic poly- 
nomials of degree £: 



H t (z,y,z) = £ C?y?{r). (B.ll) 



It follows from its definition (|B.6j) that the regular solid harmonic y™(r) 
is a homogeneous harmonic polynomial of degree I in the Cartesian compo- 
nents of r = (x,y,z) [E3 p. 71]: 



yr(r) 



,1/2 



1£ 4- 1 

-±—(£ + m)l(£-m) 
Air 



E(—x — iy) m+k (x — iy) k z e m 2k 



fc>0 



C Appendix: Gaunt Coefficients 



The so-called Gaunt coefficient [35| is the integral of the product of three 
spherical harmonics over the surface of the unit sphere in R 3 : 

(4m 3 |4m 2 |4mi) = f [Y™ 3 (Q)] * Y^{£1) Y™ 1 (fi) dQ . (C.l) 
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Gaunt coefficient can be expressed in terms of Clebsch-Gordan coefficients 
[T21 Eq. (3.192)] or 3jm symbols ^3 p. 168], respectively: 



(£3m 3 \£ 2 m2\£imi) 



(2*i + 1)(2£ 2 + 1) 



1/2 



(-1)™ 3 



4tt(24 + 1) 

(24 + 1) (2*2 + 1) (2^3 + 1 



°000 u mim 2 m3 W'*) 
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h £2 

mi m 2 



1/2 



-m 3 



(C.3) 



It follows from the orthonormality of the spherical harmonics that the 
Gaunt coefficients linearize the product of two spherical harmonics: 

^max 

Y™^n)Y™*(Q) = ^2 i2) ^m 1 + m 2 \£ 1 m 1 \£ 2 m 2 )Y e mM (n). (C.4) 



The symbol indicates that the summation proceeds in steps of two. 

The summation limits in (|C.4|) . which follow from the selection rules of the 
3jm symbols, are given by [1281 Eq. (3.1)] 



A, 



A m j n , if *max + A m j n is even , 

Amin + 1 , if *max + A min is odd , 

max(|*i - £ 2 \, \mi + m 2 \) . 



(C.5a) 
(C.5b) 
(C.5c) 



Gaunt coefficients can be computed according to either ()G2|) or ()G3|) 
via the numerous explicit expressions for Clebsch-Gordan coefficients or 3jm 
symbols described in the literature. In the case of small angular quantum 
numbers, this is a satisfactory approach, but in the case of large angular 
quantum numbers, the computation of individual Clebsch-Gordan coeffi- 
cients or 3jm symbols via their explicit expressions becomes computation- 
ally costly. In addition, a catastrophic accumulation of rounding errors can 
easily happen. A much better approach consists in the use of a recurrence 
formulas for 3jm symbols derived by Schulten and Gordon j^Tj. Although 
this recursion is neither stable upwards nor downwards, it is nevertheless 
possible to compute in this way whole strings of 3jm symbols efficiently and 
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reliably even for very large angular momentum quantum numbers |98j . In 
|128| . a FORTRAN program is described that computes simultaneously all 
Gaunt coefficients occurring on the right-hand side of ()C.4|) with the help 
of the recurrence formula and the computational algorithm of Schulten and 
Gordon EH1 - 

It may be of interest for the reader that several articles on Gaunt co- 
efficients have appeared in recent years 031 H371 11381 11391 11UUI 155] . Then, 
there are some very recent articles by Dunlad [22J EE1 EE] > but the objects he 
considered are not the usual Gaunt coefficients defined in ()0.1j) but rather 
generalizations that occur in the context of molecular multicenter integrals 
of spherical Gaussian functions. 

In this article, the following abbreviations will frequently be used: 



If the three orbital angular momentum quantum numbers l\ , £ 2 , and £ satisfy 
the summation limits (jC.5j) . then these quantities are always positive integers 
or zero. 
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